Cite this Article
Li Peng-Cheng, Chu Shih-I. Multiphoton quantum dynamics of many-electron atomic and molecular systems in intense laser fields. Chinese Physics B, 2020, 29(8): 083202  
Permissions
Multiphoton quantum dynamics of many-electron atomic and molecular systems in intense laser fields
Li Peng-Cheng1, 2, †, Chu Shih-I3
Research Center for Advanced Optics and Photoelectronics, Department of Physics, College of Science, Shantou University, Shantou 515063, China
Key Laboratory of Intelligent Manufacturing Technology of MOE, Shantou University, Shantou 515063, China
Department of Chemistry, University of Kansas, Lawrence, Kansas 66045, USA

 

† Corresponding author. E-mail: pchli@stu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11674268 and 11764038), the Natural Science Foundation of Guangdong Province, China (Grant No. 2020A1515010927), and Department of Education of Guangdong Province, China (Grant Nos. 2018KCXTD011 and 2019KTSCX037).

Abstract

We present the recent new developments of time-dependent Schrödinger equation and time-dependent density-functional theory for accurate and efficient treatment of the electronic structure and time-dependent quantum dynamics of many-electron atomic and molecular systems in intense laser fields. We extend time-dependent generalized pseudospectral (TDGPS) numerical method developed for time-dependent wave equations in multielectron systems. The TDGPS method allows us to obtain highly accurate time-dependent wave functions with the use of only a modest number of spatial grid point for complex quantum dynamical calculations. The usefulness of these procedures is illustrated by a few case studies of atomic and molecular processes of current interests in intense laser fields, including multiphoton ionization, above-threshold ionization, high-order harmonic generation, attosecond pulse generation, and quantum dynamical processes related to multielectron effects. We conclude this paper with some open questions and perspectives of multiphoton quantum dynamics of many-electron atomic and molecular systems in intense laser fields.

PACS: ;32.80.Rm;;33.80.Rv;;42.50.Hz;;42.65.Re;
Keyword:;multiphoton ionization;;above-threshold ionization;;high-order harmonic generation;;attosecond pulse generation;
1. Introduction

The fundamental phenomena from a laser beam focused on a atomic or molecular target, such as multiphoton ionization (MPI), above-threshold ionization (ATI), high-order harmonic generation (HHG), attosecond pulses generation (APG), etc., have been largely studied over the last decade. Driven by the experimental works, theoretical non-perturbative methods and various models have been proposed for the atomic and molecular dynamical behaviors in intense laser fields. The semiclassical three-step model developed by Corkum[1] and Kulander[2] has shown a fairly clear physical picture of laser-driven electronic rescattering and has been successfully extended to include some quantum effects.[3] However, with rapidly developing laser techniques, a number of very-high-order nonlinear optical processes, such as many-electron correlation effects, non-sequential double ionization, etc., have been observed from many-electron systems in intense laser pulses.[4] To study such strong nonlinear optical phenomena observed in many-electron systems, a quantitative description of the dynamics of the system will be most reliably obtained from integration of the full-dimensional time-dependent Schrödinger equation (TDSE) based on the full many-body Hamiltonian.

For the simplest many-electron system, such as helium atom with two electrons, theoretical description of this process is a six-dimensional numerical integration of the TDSE, several approaches have been proposed.[5,6] However, fully converged N-electron systems (N > 2) calculation with 3N-dimensional TDSE is still difficult to achieve and remains a major computational challenge in the study of strong-field atomic, molecular, and optical physics today. One of the approximations commonly used for the treatment of strong-field processes in the past decade is the single-active-electron (SAE) model with frozen nuclei.[7] The SAE approach is a theoretical model frequently employed to investigate scenarios in which inner-shell electrons may productively be treated as frozen spectators to a physical process of interest, providing useful insights regarding strong-field systems dynamics. However, within the SAE approach, the effects of many-electron correlation and the individual spin–orbital contribution can not be explicitly treated.

In current frontiers of laser sciences, studies of interactions between a highly-intense, ultrashort laser pulse and many-electron systems have been attracting much attention. It is desirable to explore more comprehensive formalisms for detailed treatment of strong-field processes, taking into account both the electron correlation and the structure of the excited states, and at the same time allowing for the core excitation of the many-electron systems. Density functional theory (DFT) provides a first-principles approach for the quantum-mechanical description of electrons and bypasses the need for direct calculation of a many-electron wavefunction, which is based on the earlier fundamental work of Hohenberg and Kohn[8] and Kohn and Sham.[9] It is one of the most widely used methods for ab initio calculations of the structure of atoms, molecules, and crystals.[10–14] In Kohn–Sham DFT formalism, 3N-dimensional numerical integration problem is decomposed into a set of orbitals by the electron density, leading to a set of one-electron Schrödinger-like equations to be solved self-consistently. The Kohn–Sham equations are structurally similar to the Hartree–Fock equations but include in principle exactly the many-body effects through a local exchange–correlation (xc) potential. Thus DFT is computationally much less expensive than the traditional ab initio many-electron wavefunction approaches and this accounts for its great success for large systems.[15,16]

The central theme of DFT is possible and beneficial to obtain the xc-energy functional for electronic systems, its formally exact expression is unknown. The xc-energy functional can be approximated by the local density approximation (LSDA)[17] or the generalized gradient approximation (GGA).[18–21] Such approximations are widely and successfully used to predict, understand, and design physical and chemical phenomena associated with molecules and materials. However, the xc potentials derived from these GGA energy functionals suffer similar problems as in LSDA and do not exhibit the asymptotically correct long-range Coulombic tail –1/r behavior that is expected from general considerations. The problem of the incorrect long-range behavior of the LSDA and GGA energy functionals is so-called self-interaction error (SIE).[22] Thus, whereas widely used GGA approximations allow for rather accurate predictions of the total energies of the ground states of atoms and molecules, but the excited-state energies and the ionization potentials obtained from the highest occupied orbital energies of atoms and molecules are far from satisfactory.[23] For quantitative treatment of multiphoton ionization processes, it is necessary to extend the DFT to properly account for the long-range xc potential such that both the ionization potential and the excited-state properties can be described more accurately. A Krieger–Li–Iafrate (KLI)[24,25] semianalytic treatment of the optimized effective potential (OEP)[26,27] formalism along with the use of an explicit self-interaction-correction (SIC)[22] term has been proposed independently by Chen et al.,[28] but it has been applied to the study of ground-state energies and ionization potentials only. In our previous works,[23] we have extended such a KLI-SIC technique to the time-dependent domain. In this proposed KLI-SIC procedure, similar to the original KLI method, it also allows the construction of a self-interaction-free effective potential that is orbital independent. This avoids the problems associated with the conventional SIC procedure discussed earlier in Ref. [22].

The SIC is very much needed in time-dependent density functional theory (TDDFT) for proper treatment of atomic and molecular time-dependent dynamics such as collisions or MPI processes, etc. The TDDFT is a nontrivial extension of the steady-state DFT of Hohenberg, Kohn, and Sham[8,9] to the time domain. The Runge–Gross (RG) theorem[29] provides the formal foundation of TDDFT. In RG theorem, for any interacting many-particle quantum system subject to a given time-dependent potential, all physical observables are uniquely determined by knowledge of the time-dependent density and the state of the system at any instant in time.[30] The central result of the TDDFT is a set of time-dependent Kohn–Sham equations that are structurally similar to the time-dependent Hartree–Fock (HF) equations but include in principle all the many-body effects through a local time-dependent xc potential. Telnov et al.[31] have presented an alternative nonperturbative formulation of TDDFT based on the extension of the generalized Floquet formalism,[32,33] allowing exact transformation of the time-dependent Kohn–Sham equations into an equivalent time-independent Floquet matrix eigenvalue problem. Such a TDDFT-Floquet formalism provides a time-independent approach for nonperturbative treatment of multiphoton processes of many-electron systems in the presence of intense periodic or multicolor laser fields. Tong et al.[34] have extended this time-dependent approach to the numerical solution of time-dependent Kohn–Sham-like equations in arbitrarily time-dependent fields.

In this paper, we briefly describe several new developments in TDSE and TDDFT for nonperturbative treatment of multiphoton dynamics and nonlinear optical processes of many-electron atoms and molecules in intense laser fields. In Section 2, we introduce the basic concepts and theorems of TDSE, DFT, and TDDFT. In Section 3, the multiphoton quantum dynamics of many-electron atoms and molecules are discussed, including MPI, ATI, HHG, and APG. Section 4 contains the concluding remarks. Atomic units are used in this paper.

2. Theoretical methods
2.1. Time-dependent Schrödinger equation

Consider an N-electron atom (ion) in an electromagnetic field described by the Coulomb gauge obeys the TDSE, if we neglect relativistic effects and nuclear effects when the atoms interact with strong laser fields, the TDSE can be written in the dipole approximation as (in atomic units)[35]

where X ≡(q1,q2,…,qN) denotes the ensemble of the coordinates of the N electrons and qi ≡ (ri,σi) are the space and spin coordinates of the i-th electron, A(t) is the vector potential, and

is the total momentum operator. H0 is the time-independent Hamiltonian of the N-electron atom (ion) in the absence of the electromagnetic field, and given by

where V denotes the sum of all the interactions within the atomic system in the absence of the radiation field. However, the property of gauge invariance allows us to simplify the TDSE by an appropriate choice of gauge. In the velocity gauge (VG), the TDSE can be be written as[35]

In the length gauge (LG), the TDSE can be be written as

where E(t) is the electric field and

is the sum of the coordinates of the N electrons. The formal theoretical identity of results obtained in the two gauges does not necessarily mean that such an identity can be obtained in practical numerical calculations. Since quantum mechanics is gauge invariant, the results obtained from the two gauges after the end of the induced electron motion should be the same, but two formulations may impose different demands on computations in terms of basis size and cost of CPU time, the previous study[36] demonstrates that the propagation of TDSE solved in the LG for nonperturbative treatment of atoms and molecules in intense laser fields is faster than the case in the velocity gauge, and the VG is the optimal gauge for the ATI calculation. Consider hydrogenic atoms and ions in intense laser fields, corresponding to the case N = 1, we can respectively rewrite the TDSE in the LG and the VG in the forms

where

In recent years, advances in computer technology have allowed the numerical solution of the TDSE for two-electron atoms to be obtained, but the generalization of this approach to the solution of the TDSE for multielectron atoms is beyond the present computer capabilities.[5,6] One of the approximations commonly used for the treatment of strong-field processes of the complex atom is the SAE model with frozen core.[7] The SAE approach has been successfully applied to the investigation of HHG of atoms and molecules in intense laser fields, providing useful insights regarding strong-field atomic and molecular dynamics. In the SAE model, the time-independent Hamiltonian H0 can be written in the form

where V(r) is the atomic model potential. The SAE model is capable of describing the physics of multiphoton single ionization and more generally multiphoton processes involving single one-electron transitions. For the complex atom with a nucleus of atomic number Z and N electrons, a standard approximation used to perform time-independent calculations for many-electron atomic systems is the time-dependent Hartree–Fock (TDHF) approximation method, based on the independent particle model, in which each of the atomic electrons moves in a self-consistent field that takes into account the attraction of the nucleus and the average effect of the repulsive interactions due to the other electrons.

2.2. Time-dependent density function theory

The TDDFT is an extension of the DFT[8,9] to the time domain, providing an alternative approach for nonperturbative treatment of multiphoton processes of many-electron systems in intense laser fields. The central theme of modern TDDFT is a set of time-dependent Kohn–Sham (TDKS) equations which are structurally similar to the TDHF equations but include in principle all the many-body effects through a local time-dependent xc potential. Most applications of TDDFT[37,38] fall in the regime of linear or nonlinear response in weak fields for which the perturbation theory is applicable. In our previous works,[34] we have proposed a TDDFT with optimized effective potential (OEP) and SIC developed for nonperturbative treatment of intense-field atomic multiphoton processes of the atomic and molecular systems. The TDDFT/OEP-SIC formalism takes into account both the electron correlation effect and the structure of the excited states, and at the same time allows for the core excitation of the many-electron systems in intense laser pulses. The method takes into account the dynamic response of all the electronic shells to the external fields and has been applied successfully to the nonperturbative study of MPI and HHG of atoms and diatomic molecules[39,40] in intense laser fields. In the TDDFT/OEP-SIC frame, the TDKS equations of N-electron systems in intense laser fields can be written as

where is the time-dependent OEP with SIC depending upon the total electron density ρ(r,t), Nσ ( = N↑ or N↓) is the total number of electrons for a given spin σ, and the total number of electrons in the system is N = ∑σNσ. The total electron density ρ(r,t) is determined by the single-electron orbital wave functions ψi\sigma(r,t) as

The time-dependent effective potential can be written in the following general form:

where

is the electron–electron Coulomb interaction energy, υext(r,t) is the external potential due to the interaction of the electron with the external laser field and the nucleus, and is the time-dependent exchange–correlation potential with SIC. It has the following form:

where

2.3. Generalized pseudospectral method

The generalized pseudospectral (GPS) method can be applied to the calculation of multiphoton dynamics and nonlinear optical processes of many-electron atoms and molecules in intense laser fields.[41,42] The essence of the GPS method is to map the semi-infinite domain [0, ∞] or [0, rmax] into the finite domain [–1,1] using a non-linear mapping r = r(x), followed by the Legendre or Chebyshev pseudospectral discretization.[43] This allows for denser grids near the origin, leading to more accurate eigenvalues and eigenfunctions and the use of a considerably smaller number of grid points than those of the equal-spacing grid methods. For example, in Table 1, we show the bound state energies of calculated by the GPS method,[44] only 72 × 24 grid points for the spatial coordinates were used, we were able to obtain the energies of the first several bound states of with very high accuracy. Thus the GPS method with the parameters given above can reproduce the initial states for the time propagation procedure as well as the propagation matrices with sufficiently high accuracy.

Table 1.

High-precision bound state energies of at internuclear separation 2 a.u.

.

To the time propagation of the wave function, we have extended the GPS method to the time-dependent generalized pseudospectral (TDGPS).[45] The numerical scheme of the TDGPS method consists of two essential steps: (i) The spatial coordinates are optimally discretized in a nonuniform fashion by means of the denser grids near the nuclear origin and sparser grids for larger distances. (ii) A second-order split-operator technique in the energy representation, which allows the explicit elimination of undesirable fast-oscillating high-energy components, is used for the efficient and accurate time propagation of the wave function. To introduce the numerical process of TDGPS in detail, we consider a hydrogen atom case which is convenient due to the similar procedure for many-electron atomic and molecular systems in intense laser fields. The solution of TDSE in Eq. (7) can be written as

where are the spherical harmonics and are the associated Legendre polynomials. The TDSE becomes

Based on the GPS numerical technique, we employ the following coordinate transformation r ∈ [0,rmax ] to x ∈ [ –1,1], a suitable algebraic mapping for atomic structure calculations is provided by the following form:

where L and α = 2L/rmax are the mapping parameters. The radial function in Eq. (21) is interpolated at the Legendre–Lobatto collocation point {x},

where is the first derivative of the N-th order Legendre polynomial PN, and its quadrature weight is

Let us introduce the following discretized wavefunction:

where

We have extended the second-order split-operator technique in spherical coordinates[47–49] for the time propagation of the TDSE

The stationary part of the Hamiltonian is operated as

where the evolution matrix Sij(l) is constructed from the eigenstates of the Hamiltonian,

such that

Once the time-dependent wave function Ψ(r,t) is available by solving the TDSE or TDDFT/OEP-SIC equations, it allows us to observe the multiphoton quantum dynamics of many-electron atomic and molecular systems in intense laser fields.

3. Multiphoton quantum dynamics of atoms and molecules
3.1. MPI

Traditionally many theoretical studies of MPI, ATI, and HHG processes in many-electron atoms and molecules are based on the strong-field approximation (SFA),[3] this approach has its origin in earlier works of Keldysh, Faisal, and Reiss[50–52] as well as in the semiclassical rescattering model.[1,2] While SFA-based models result in rather simple theoretical expressions and produce electron spectra that qualitatively resemble those obtained with accurate numerical wave functions, they fail to give quantitative agreement with more accurate theories. Many attempts have been made recently to improve SFA,[53–55] it has intrinsic restrictions and cannot compete with ab initio calculations for accuracy of the results. On the other hand, SFA-based theories neglect the multielectron dynamics of the target systems. However, the multielectron effects due to the electron exchange and correlation may be significant even when the inner electrons are strongly bound and are not excited by the driving laser field.[56,57] In our previous work,[34,39,40,58] we have extended the TDDFT with proper long-range potential to an all-electron three-dimensional (3D) ab initio study of the MPI of atoms and molecules in intense laser fields, a subject of much current experimental interests.[59–61] Consider a many-electron Ar atom, the time-dependent wave function ψiσ(r,t) in Eq. (11) is obtained by the TDGPS method, so the total electron density ρ(r,t) can also be determined. The time-dependent muitiphoton ionization probability of the i-th spin orbital can be calculated according to

where

is the time-dependent population survival probability of the i-th spin-orbital. The molecular muitiphoton ionization probability calculated by the TDDFT has a similar definition.[39]

In Fig. 1(a) we present the time-dependent ionization probabilities of the Ar atom from individual valence electron orbital 3s, 3p0, and 3p1 calculated by solving the TDDFT/OEP-SIC equations. In calculation, the laser wavelength is 800 nm with the laser peak intensity I = 8.0 × 1013 W/cm2, and the pulse length is 20 optical cycles with sin2 pulse shape. It is clear that the ionization probability of 3p0 is the maximum one, and the 3p1 subshell ionization rate is higher than that of the 3s subshell. According to our previous works,[62] the np0 valence electron will always be the easiest one to ionize, while the ionization probability of the ns and np1 electrons will depend upon the relative importance of the binding energy and orientation factors, as well as the possible enhancement due to accidental multiphoton resonant excitation. Since the laser peak intensity plays an important role in the ionization behavior of the valence electrons and the generation of harmonic spectrum, we present the time-dependent ionization probabilities of the Ar atom from individual valence electron subshells 3s, 3p0, and 3p1 by driving laser pulse with the laser peak intensity I = 3.0×1014 W/cm2, as shown in Fig. 1(b), other laser parameters used are the same as those in Fig. 1(a). The results show that the time-dependent ionization probabilities of the Ar atom from individual valence electron subshells 3s, 3p0, and 3p1 are increased, and it is clear that the total time-dependent ionization probability is larger than 20%. The results indicate that the time-dependent ionization probability plays an important role in the propagation dynamics of the laser pulse.

Fig. 1. (a) Time-dependent ionization probability of 3s, 3p0, and 3p1 valence electrons of Ar atom in intense ultrashort laser fields with the laser peak intensity I = 8.0 × 1013 W/cm2. The laser wavelength is 800 nm, and the pulse length is 20 optical cycles with sin2 pulse shape. (b) Same as (a) for the laser peak intensity I = 3.0×1014 W/cm2.

In Figs. 2(a) and 2(b), we present the time-dependent ionization probability of the different molecular orbital electrons of N2 and CO in 800 nm, sin2 laser with 20 optical cycles in pulse duration, and the peak intensities of laser fields used are I = 8×1013 W/cm2 and I = 5.0×1013 W/cm2, respectively. The laser field is assumed to be parallel to the internuclear axis, and the internuclear distance for the N2 (R = 2.072 a.u.) and CO (R = 2.132 a.u.) molecules is fixed at its equilibrium distance. The orbital structure and ionization potentials of the two molecules under consideration are close to each other. That is why one can expect similar behavior in the laser field with the same wavelength. The ground state electronic configuration of N2 and CO is and , respectively. The multiphoton ionization in the laser field is dominated by HOMO, that is in N2 and 5σ in CO. The smaller the ionization potential of the electronic shell is, the easier it can be ionized. That is why HOMO is generally expected to give the main contribution to the MPI probability. However, the orbital probabilities of HOMO-1 (1π) of N2 and CO have different behaviors with the laser intensity. The reason is that the N2 molecule is symmetric with respect to inversion. On the contrary, the CO molecule has a permanent dipole moment, the multiphoton ionization depends on the direction of the external field with respect to the position of the carbon and oxygen nuclei.

Fig. 2. Time-dependent ionization probability of the different molecular orbitals electrons of N2 and CO in 800 nm, sin2 laser pulse with 20 optical cycles in pulse duration. (a) N2 molecule, for the peak intensity of laser field I = 8×1013 W/cm2. (b) CO molecule, for the laser peak intensity I = 5.0×1013 W/cm2.

In Figs. 3(a) and 3(b), we present the orientation-dependent MPI probabilities for N2 molecule and Ar atom at the peak intensities of 1×1014 W/cm2 and 5×1014 W/cm2, respectively. We used the laser wavelength 800 nm and the sine-squared envelope with 20 optical cycles. The orientation dependence of the total MPI probability is in a good accord with the experimental observations.[59,60] The maximum MPI probability for N2 molecule corresponds to the parallel orientation and reflects the symmetry of the HOMO of N2. However, multielectron effects are quite important for N2, particularly at intermediate orientation angles.[39] Despite the orbital probabilities have local minima and maxima, the total probability shows monotonous dependence on the orientation angle. With increasing the peak intensity of the laser field, the orientation angle distribution of the total ionization probability becomes more isotropic. For comparison, we also show the ionization probability of the Ar atom, whose ionization potential is similar to that of N2 (HOMO). As one can see from Figs. 3(a) and 3(b), the absolute values of the ionization probabilities of N2 and Ar are close to each other. However, the inner shell contributions are less important for Ar, the total probability is dominated by the highest-occupied (3p) shell contribution.

Fig. 3. The orientation-dependent MPI probabilities for N2 molecule and Ar atom at the peak intensity (a) I = 1×1014 W/cm2 and (b) I = 5×1014 W/cm2.
3.2. ATI

The investigation of above-threshold ionization (ATI)[63] has benefited from the rapid progress of the laser technology. Grasbon et al.[64] have measured ATI photoelectron spectra for noble gas atoms ionized with intense few-cycle laser pulses. Resolutions of the photoelectron momentum distribution in experiments are unprecedentedly high, which has made detailed comparison with theoretical calculations possible. Since the orbital-dependent measurements of the p-state photoelectron momentum distribution are already available,[65] theoretical calculations beyond the SAE approximation are urgently needed.

We have extended the TDDFT/OEP-SIC for the calculation of the photoelectron momentum distribution of the noble gas atoms driven by linearly polarized laser fields. For the ATI calculation, the wave function ψiσ(r,t) is split into inner and outer regions by a smooth masking function, and the photoelectron momentum distribution is found from the outer-region wave function that is propagating in the momentum space with the Volkov Hamiltonian in the velocity gauge.[66–68] We study the photoelectron-momentum-distribution of each individual electron at the end of the time evolution t = tf, given by

where is the Fourier transform of the outer-region wave function, as well as the photoelectron-momentum-distribution of all electrons, which can be written as

Figures 4(a) and 4(b) show the photoelectron-momentum-distributions of the 2p and 3p states of Ne and Ar, respectively. In calculation, the noble gas atoms are driven by the 800-nm, linearly polarized, cos2 pulse envelope, 20-cycle laser pulse along the z axis with a peak intensity of I = 1.0×1014 W/cm2. The results indicate that the low-energy photoelectron-momentum-distributions of the p-state electrons of Ne and Ar are different along the z axis, and the bound states and subshell structures play an important role in the photoelectron-momentum distributions.

Fig. 4. The p-state photoelectron-momentum-distribution cross section of (a) Ne and (b) Ar driven by the 800-nm, linearly polarized, 20-cycle laser pulse along the z axis with a peak intensity of I = 1.0×1014 W/cm2.
3.3. HHG
3.3.1. Resonant enhancement

Consider a helium atom, to obtain the accurate calculation of the harmonic spectra of He, an angular-momentum-dependent model potential is constructed in the following form:[69–71]

where α is the He+ core dipole polarizability, W6 is a core cutoff function[72] given by

and rc is an effective He+ core radius.

We assume that the laser polarization is along the z-axis. The TDSE in Eq. (21) becomes

where V(r,t) is the time-dependent atom–field interaction, and H0 represents the unperturbed atom Hamiltonian, which is

where is the spherical harmonic. In the present work we find that it is sufficient to use five different angular-momentum-dependent model potentials, namely, all for states with l = 0,1,2,3 and l ≥ 3. The values of the parameters determined are listed in Table 2. Table 3 presents a comparison of the bound-state energies predicted by this model potential and the experimental values. The agreement is well satisfactory in all cases.

Table 2.

Model potential parameters for He (in a.u.).

.
Table 3.

Comparison of the calculated He atomic energies with the experimental values (in a.u.). For each angular momentum l, two rows of energies En,l are listed: the upper row refers to the calculated model-potential energies, and the lower row refers the experimental values.[NIST]

.

The TDSE is solved accurately and efficiently by means of the TDGPS. Once the time-dependent wave function ψ(r,t) is available, we can calculate the expectation value of the induced dipole moment,

From Ehrenfest theorem, the single-atom harmonic spectra can be calculated by the Fourier transformation of the time-dependent dipole moment. We employ the widely-used semi-classical approach, where the basic expressions come from classical electrodynamics. The spectral density of the radiation energy emitted is given by the following expression:[74]

where c is the speed of light, and is the Fourier transform of the time-dependent dipole acceleration d(t).

In Figs. 5(a) and 5(b), the HHG of He atom driven by a few-cycle 760-nm laser pulse as a function of the laser intensity and the photon energy in the single atom response are shown. To simulate the experimental condition reported in Ref. [75], we use the sine-squared envelope, and 4 optical cycles (FWHM pulse duration is ∼ 5 fs). The main pulse sits on a pedestal of a weaker field (4% of the peak intensity) which also has a sine-squared shape and duration of 20 optical cycles. Thus the stronger field with 4 optical cycles is preceded by 8 optical cycles and also followed by 8 optical cycles of the weaker field. The result shows that the yields of the resonance-enhanced harmonic spectra strongly depend on the laser peak intensity, and the yield of the 13th harmonic with the 760-nm laser frequency, which is just on resonant with the transition of the 1s2-1s2p1P (21.22 eV), is always largely enhanced, and the yields of the 15th harmonic only have a slight enhancement with the laser intensity. Furthermore, a maximum yield of HHG is located at I = 9.0×1013 W/cm2.

Fig. 5. (a) Resonance-enhanced structures of the HHG from He as a function of the laser intensity and the photon energy in the single atom response. (b) Same as (a) but for I = 7.0×1013 W/cm2, I = 8.0×1013 W/cm2, and I = 9.0×1013 W/cm2, respectively. The black dashed lines in (a) and (b) indicate the transition energy of the 1s2-1s2p1P (21.22 eV), and the black solid lines indicate ionization potential Ip.

The He atom only has two electrons, all occupying s states, the resonance-enhanced HHG described by the TDSE with an angular-momentum-dependent model potential is sufficient for the single-atom response. For the complex atoms and molecules, the TDDFT can obtain more accurate calculation of the harmonic spectra. In the following, we extend the TDGPS procedure to the numerical solution of the TDKS equations for N-electron atomic systems in intense laser fields. Once the time-dependent wave function ψiσ(r,t) in Eq. (11) is obtained, the total electron density ρ(r,t) can also be determined. The induced dipole moment can be expressed as

where niσ is the electron occupation number. However, the laser field interacts with the atoms or molecules, a full description of experimentally observed HHG spectra requires not only the theoretical treatment of the nonlinear laser–atom interaction but also the macroscopic propagation of the radiation through the nonlinear optical medium. The propagation of the laser-field EL and harmonic-field EH in a macroscopic medium is described by 3D Maxwell wave equation, it can be written in the following form:[76–80]

where Pl(r,t) is the linear polarization, which has the following form:

here the linear susceptibility χ includes the dispersion δ and absorption effects β in the medium. The dispersion δ can be obtain from Ref. [17] and particularly it is respectable for the main contributions of macroscopic resonance-enhanced HHG. The absorption effect is not important in phase-matched below-threshold harmonic generation, so it is neglected. The nonlinear polarization component Pnl(r,t) can be written as

where n0 is the neutral atom density, ne(r,t) is the free-electron density, and d(r,t) is the single-atom induced dipole moment calculated according to Eq. (11). Equations (42) and (43) are solved in the frequency domain.

In fact, when the pressure and the laser intensity are low enough, the incident laser field is not modified in the medium, and only the harmonic field has to be propagated. Once equation (43) is solved, the macroscopic HHG spectrum can be obtained by integration of the signal in the plane perpendicular to the propagation direction

Since the TDDFT calculation for single-atom response for feeding the macroscopic propagation equations is very time consuming, the fast and efficient computing is desirable. For this purpose, we design a parallel algorithm based on the graphics processing unit (GPU) machines for the calculation of the TDDFT plus Maxwell wave equation.

In Fig. 6(a), we present the phase-matching of resonance-enhanced structures (RESs) from Ar atom simulated by the TDDFT/OEP-SIC including the propagation effects. The laser peak intensity in the center of gas jet is I = 3.0×1013 W/cm2, the confocal parameter of the laser beam is 35 mm, the length of gas jet is 2.0 mm, the pressure is 40 Torr, and the target is set at the laser focus –2 mm. The harmonic spectrum is characterized by coherent line emissions which are located at the energy range of the labeled 3p6 to 3p5ns and 3p5nd resonance states near the 9th harmonic. The harmonic emission in argon indicates that the phase-matching process results in spectrally narrow, energy-shifted RESs when compared with the single-atom emission. The result shows that the phase-matching occurs primarily in the vicinity of the resonances, in good agreement with the experiment.[75] Figure 6(b) shows the phase-matching of resonance-enhanced structures from Ar atom for the 60 Torr pressure case. As one can see that the resonance-enhanced structure appears more fine peaks near the 9th harmonic order.

Fig. 6. (a) Phase-matching of resonance-enhanced structures from Ar atom simulated by the TDDFT/OEP-SIC in a 730-nm laser field. The shift of the phase-matched RESs relative to the bound state energies is indicated in the inset. The laser peak intensity in the center of gas jet is I = 3.0×1013 W/cm2, the confocal parameter of the laser beam is 35 mm, the length of gas jet is 2.0 mm, the pressure is 40 Torr, and the target is set at the laser focus –2 mm. (b) Same as (a) for the 60 Torr pressure.
3.3.2. Phase-matching of HHG

Phase matching plays an essential role in determining the efficiency of macroscopic HHG. The condition for perfect phase-matching (PM) in HHG is obtained when the wavevector of the q-th harmonic kq equals the dipole wavevector kd, which can be written as[81]

where ϕ(r,z,t) and Φ(r,z,t) are the phase accumulated during propagation and the phase of the driving field, respectively. Both the phases will finally determine how the HHG emission coherently builds up so that, in general, in order to achieve efficient HHG along a defined propagation direction phase, matching requirements must be fulfilled.[82,83]

As is well known, such a condition is fulfilled in a Gaussian beam, on-axis and after the focus, when both the phase and intensity decrease along the propagation direction. When the gas medium is placed before the focus, the PM condition is only fulfilled in restricted regions off-axis where the phase and intensity gradients in the radial direction become important. In Figs. 7(a) and 7(c), we show the evolution of harmonic intensities in space for the 7th order harmonic (7H) and the 15th order harmonic (15H) of H2 molecules obtained from the mac TDDFT in an 800-nm laser field, the corresponding phase difference of the HHG is presented in Figs. 7(b) and 7(d), respectively. The laser peak intensity in the center of gas jet is I = 8.0×1013 W/cm2, the pulse length is 20 optical cycles with cos2 pulse shape, the confocal parameter of the laser beam is 25 mm, the length of gas jet is 4.0 mm, the pressure is 10 Torr, and the target is set at the laser focus. It is clear that the 7H and 15H can manipulate the gas jet position to obtain optimum harmonic yields, and the phase-matching condition is fulfilled in restricted phase difference from –π/2 to π/2 (see the green solid lines).

Fig. 7. (a) and (c) Evolution of harmonic intensities in space for the 7th order harmonic (7H) and the 15th order harmonic (15H) of H2 molecules obtained from the mac TDDFT in an 800-nm laser field. The laser peak intensity in the center of gas jet is I = 8.0×1013 W/cm2, the pulse length is 20 optical cycles with cos2 pulse shape. (b) and (d) The corresponding phase difference of the HHG.

Recently, Willner et al.[84] have developed a novel dual-gas quasiphase matching (QPM) concept based on alternating a HHG generating medium with passive matching hydrogen zones. The choice condition of the dual-gas target proposed by their previous results is that the HHG medium must have a higher ionization potential than the phase matching medium. In Fig. 8, we show a scheme of the dual-gas QPM in HHG, but the helium atom which only acts as the phase-matching medium has a larger ionization potential than those argon HHG medium. In Fig. 9, the enhancement of the HHG calculated by the mac TDDFT is presented due to the dual-gas QPM. For reference, the corresponding one-gas result (blue dotted lines) is also presented. In calculation, the laser wavelength is 800 nm, the laser peak intensity in the center of gas jet is I = 2.3×1014 W/cm2. For one-gas jet, each gas jet d1 = 0.1 mm (Ar), the total gas length is 0.6 mm, and the pressure is 20 Torr. For two-gas jet, each gas jet d1 = 0.1 mm (Ar) and d2 = 0.1 mm (He), the total gas length is 0.6 mm for Ar and 0.5 mm for He. The pressure is 20 Torr for Ar and 80 Torr for He. The gas get is put 2 mm after the focus. The results show that the HHG in dual-gas QPM is largely enhanced.

Fig. 8. A scheme of dual-gas QPM
Fig. 9. Enhancement of the HHG in dual-gas QPM. For the reference, the corresponding one-gas result (blue dotted lines) is presented.
3.4. APG

The study of the generation of the APG is a subject of much interest in ultrafast science and technology in the last decade.[85] Currently, the production of the APG by means of the superposition of broadband supercontinuum in HHG is one of the most promising routes. For the generation of the broadband supercontinuum harmonic spectra and ultrashort isolated attosecond pulse, several techniques have been proposed, such as the availability of few-cycle laser pulse,[86,87] double optical gating,[88] multi-cycle driver laser pulses directly from an amplifier,[89] and two- or multi-color laser scheme.[90–99] As an example, consider a helium atoms driven by two color combined laser field, the HHG spectrum can be obtained by the Fourier transformation of time-dependent dipole moment d(t), given by

By superposing several harmonics, an ultrashort pulse can be generated

Here, we choose two color combined laser fields in the follow form:

where E1 and E2 are the amplitudes with the Gaussian durations f1(t) and f2(t), respectively, and ω1 and ω2 are the corresponding frequencies. ϕ is the carrier envelope phase, and td is the time delay between the two laser pulses.

Figure 10(a) shows the HHG spectrum from He atoms driven by the two-color combined laser field. In our calculation, ω1 = 2400 nm, ω2 = 1200 nm, E1 = 0.119 a.u., E2 = 0.053 a.u., and td = –0.38T1 (T1 = 2π/ρ1). It is evident that there is a two-plateau structure with the cutoffs at the 800th order and 2000th order, respectively, and the second plateau is smooth, which is advantageous to produce the single ultrashort attosecond. In Fig. 10(b), we show an isolated 26 as pulse by superposing the harmonics from the 1500th to the 2000th order. This attosecond pulse is regular, and the duration is close to one atomic unit of time.

Fig. 10. (a) High-order harmonic spectra of He atoms driven by two color combined laser field with the proper time delay td. (b) Attosecond pulse generation by superposing hundreds of the harmonics shown in (a).
4. Conclusion and outlook

In summary, we have presented TDSE and TDDFT approaches recently developed for accurate and efficient treatment of the time-dependent dynamics of many-electron atomic and molecular systems. Both the TDSE and TDDFT can be solved by means of the TDGPS methods. The generalized pseudospectral technique allows the construction of non-uniform and optimal spatial grids, denser mesh nearby each nucleus and sparser mesh at longer range, leading to high-precision solution of both electronic structure and time-dependent quantum dynamics with the use of only a modest number of spatial grid points. The TDSE and TDDFT formalism along with the use of the time-dependent GPS numerical technique provides a powerful new nonperturbative time-dependent approach for exploration of the electron correlation and multiple orbital effects on strong field processes. The procedure is demonstrated by several case studies of MPI, ATI, HHG, and APG of atomic and molecular systems, such as He, Ne, and Ar atoms, H2, N2, and CO molecules. The TDDFT is the primary approach available for the treatment of time-dependent processes of many-electron quantum systems in strong fields. Further extension of the self-interaction-free TDDFT approaches to larger molecular systems will be valuable and can lead to significant advancement in the understanding of strong-field chemical physics and atomic and molecular physics in the future.

Reference
[1] Corkum P B 1993 Phys. Rev. Lett. 71 19941
[2] Kulander K C Schafer K J Krause J L 1993 Proceedings of the Workshop on Super-Intense Laser Atom Physics (SILAP) III New York Plenum Press 316 95 https://link.springer.com
[3] Lewenstein M Balcou Ph Ivanov M Yu L’Huillier A Corkum P B 1994 Phys. Rev. 49 2117
[4] Popmintchev T et al. 2012 Science. 336 1287
[5] Parker J S Moore L R Dundas D Taylor K T 2000 J. Phys. B: At. Mol. Opt. Phys. 33 L691
[6] Nikolopoulos L A A Lambropoulos P 2007 J. Phys. B: At. Mol. Opt. Phys. 40 1347
[7] Krause J L Schafer K J Kulander K C 1992 Phys. Rev. 45 4998
[8] Hohenberg P Kohn W 1964 Phys. Rev. 136 B864
[9] Kohn W Sham L J 1965 Phys. Rev. 140 A1133
[10] Parr R Yang W 1989 Density-Function Theory of Atoms and Moleculesm New York Oxford University Press 10.1002/qua.560470107
[11] Dreizler R M Gross E K U 1990 Density Functional Theory, An Approach to the Quantum Many-Body Problem Berlin Springer-Verlag 10.1007/978-3-642-86105-5
[12] Gross E K U Dreizler R M 1995 Density Functional Theory, NATO Advanced Study Institute, Series B: Physics New York Plenum
[13] March N H 1992 Electron Density Theory of Atoms and Molecules San Diego Academic 10.1021/j100209a022
[14] Labanowski J K Andzelm J W 1991 Density Functional Methods in Chemistry Berlin Springer-Verlag
[15] Son S K Chu S I 2009 Phys. Rev. 80 011403R
[16] Li P C Sheu Y L Chu S I 2020 Phys. Rev. 101 011401R
[17] Vosko S H Wilk L Nusair M 1980 Can. J. Phys. 58 1200
[18] Becke A D 1988 Phys. Rev. 38 3098
[19] Lee C Yang W Parr R G 1988 Phys. Rev. 37 785
[20] Perdew J P Wang Y 1986 Phys. Rev. 33 8800
[21] Zhao Q Parr R G 1992 Phys. Rev. 46 R5320
[22] Perdew J P Zunger A 1981 Phys. Rev. 23 5048
[23] Tong X M Chu S I 1997 Phys. Rev. 55 3406
[24] Krieger J B Li Y Iafrate G J 1990 Phys. Lett. 146 256
[25] Krieger J B Li Y Iafrate G J 1992 Phys. Rev. 46 5453
[26] Sharp R T Horton G K 1953 Phys. Rev. 90 317
[27] Talman J D Shadwick W F 1976 Phys. Rev. 14 36
[28] Chen J Krieger J B Li Y Iafrate G J 1996 Phys. Rev. 54 3939
[29] Runge E Gross E K U 1984 Phys. Rev. Lett. 52 997
[30] Gross E K U Kohn W 1985 Phys. Rev. Lett. 55 2850
[31] Telnov D Chu S I 1997 Chem. Phys. Lett. 264 466
[32] Chu S I 1985 Adv. At. Mol. Phys. 21 197
[33] Chu S I 1989 Adv. Chem. Phys. 73 739
[34] Tong X M Chu S I 1998 Phys. Rev. 57 452
[35] Joachain C J Kylstra N J Potvliege R M 2012 Atoms in intense laser fields Cambridge Cambridge University Press
[36] Han Y C Madsen L B 2010 Phys. Rev. 81 063430
[37] Zangwill A Soven P 1980 Phys. Rev. 21 1561
[38] Mahan G D Subbaswamy K R 1990 Local Density Theory of Polarizability New York Plenum Press 10.1007/978-1-4899-2486-5
[39] Telnov D Chu S I 2009 Phys. Rev. 80 043412
[40] Heslar J Telnov D Chu S I 2011 Phys. Rev. 83 043414
[41] Wang J Chu S I Laughlin C 1994 Phys. Rev. 50 3208
[42] Yao G Chu S I 1993 Chem. Phys. Lett. 204 381
[43] Canuto C Hussaini M Y Quarteroni A Zang T A 1997 Spectral Methods in Fluid Dynamics Berlin Springer 10.1007/978-3-642-84108-8
[44] Telnov D Chu S I 2011 Comput. Phys. Commun. 182 18
[45] Tong X M Chu S I 1997 Chem. Phys. 217 119
[46] Madsen M M Peek J M 1971 At. Data 2 171
[47] Feit M D Fleck J A Jr Steiger A 1982 J. Comput. Phys. 47 412
[48] Jiang T F Chu S I 1992 Phys. Rev. 46 7322
[49] Hermann M R Fleck J A Jr 1988 Phys. Rev. 38 6000
[50] Keldysh L V 1964 Zh. Eksp Teor. Fiz. 47 1945
[51] Faisal F H M 1973 J. Phys. B 6 L89
[52] Reiss H R 1980 Phys. Rev. A 22 1786
[53] Madsen C B Madsen L B 2006 Phys. Rev. 74 023403
[54] Chen Z Morishita T Le A T Lin C D 2007 Phys. Rev. 76 043402
[55] Odžak S Milošević D B 2009 Phys. Rev. 79 023414
[56] Patchkovskii S Zhao Z Brabec T Villeneuve D M 2006 Phys. Rev. Lett. 97 123003
[57] Santra R Gordon A 2006 Phys. Rev. Lett. 96 073906
[58] Li P C Chu S I 2013 Phys. Rev. 88 053415
[59] Pavičić D Lee K F Rayner D M Corkum P B Villeneuve D M 2007 Phys. Rev. Lett. 98 243001
[60] Thomann I Lock R Sharma V Gagnon E Pratt S T Kapteyn H C Murnane M M Li W 2008 J. Phys. Chem. 112 9382
[61] Meckel M Comtois D Zeidler D Staudte A Pavičić D Bandulet H C Pépin H Kieffer J C Dörner R Villeneuve D M Corkum P B 2008 Science 320 1478
[62] Tong X M Chu S I 2001 Phys. Rev. 64 013417
[63] Zhou X C Shi S Li F Yang Y J Chen J Meng Q T Wang B B 2019 Chin. Phys. B 28 103201
[64] Grasbon F Paulus G G Walther H et al. 2003 Phys. Rev. Lett. 91 173003
[65] Fleischer A Wörner H J Arissian L Liu L R Meckel M Rippert A Dörner R Villeneuve D M Corkum P B Staudte A 2011 Phys. Rev. Lett. 107 113003
[66] Murakami M Zhang G P Chu S I 2017 Phys. Rev. 95 053419
[67] Tong X M Hino K Toshima N 2006 Phys. Rev. 74 031405(R)
[68] Arbó D G Yoshida S Persson E Dimitriou K I Burgdorfer J 2006 Phys. Rev. Lett. 96 143003
[69] Li P C Laughlin C Chu S I 2014 Phys. Rev. 89 023431
[70] Li P C Sheu Y L Laughlin C Chu S I 2014 Phys. Rev. 90 041401R
[71] Li P C Sheu Y L Laughlin C Chu S I 2015 Nat. Commun. 6 7178
[72] Chu X Chu S I Laughlin C 2001 Phys. Rev. 64 013406
[73] Kramida A Ralchenko Yu Reader J NIST ASD Team 2013 NIST Atomic Spectra Database (ver. 5.1) National Institute of Standards and Technology Gaithersburg, MD physics.nist.gov/asd
[74] Landau L D Lifshitz E M 1975 The classical theory of fields New York Pergamom Press
[75] Chini M et al. 2014 Nat. Photon. 8 437
[76] Gaarde M B Tate J L Schafer K J 2008 J. Phys. 41 132001
[77] Tosa V Kim H T Kim I J Nam C H 2005 Phys. Rev. 71 063807
[78] Jin C Le A T Trallero-Herrero C A Lin C D 2011 Phys. Rev. 84 043411
[79] Geissler M Tempea G Scrinzi A Schnürer M Krausz F Brabec T 1999 Phys. Rev. Lett. 83 2930
[80] Pan Y Guo F M Yang Y J Ding D J 2019 Chin. Phys. B 28 113201
[81] Vozzi C Negro M Calegari F Stagira S Kovács K Tosa V 2011 New J. Phys. 13 073003
[82] Balcou P et al. 1997 Phys. Rev. 55 3204
[83] Pfeifer T et al. 2006 Rep. Prog. Phys. 69 443
[84] Willner A et al. 2011 Phys. Rev. Lett. 107 175002
[85] Hernández-García C Pérez-Hernández J A Popmintchev T Murnane M M Kapteyn H C Jaron-Becker A Becker A Plaja L 2013 Phys. Rev. Lett. 111 033002
[86] Goulielmakis E et al. 2008 Science 320 1614
[87] Sansone G et al. 2006 Science 314 443
[88] Mashiko H Gilbertson S Li C Khan S D Shakya M M Moon E Chang Z H 2008 Phys. Rev. Lett. 100 103906
[89] Zeng B Chu W Li G H Yao J P Ni J L Zhang H S Cheng Y Xu Z Z Wu Y Chang Z H 2012 Phys. Rev. 85 033839
[90] Takahashi E J Lan P F Muecke O D Nabekawa Y Midorikawa K 2010 Phys. Rev. Lett. 104 233901
[91] Zeng Z Cheng Y Song X Li R Xu Z 2007 Phys. Rev. Lett. 98 203901
[92] Wang Z Li Y Wang S Y Hong W Y Zhang Q B Lu P X 2013 Phys. Rev. 87 033822
[93] Mauritsson J Johnsson P Gustafsson E L’Huillier A Schafer K J Gaarde M B 2006 Phys. Rev. Lett. 97 013001
[94] Lan P F Lu P X Cao W Li Y H Wang X L 2007 Phys. Rev. A 76 011402R
[95] Lan P F Takahashi E J Midorikawa K 2010 Phys. Rev. 82 053413
[96] Wu J Zhang G T Xia C L Liu X S 2010 Phys. Rev. 82 013411
[97] Li P C Liu I L Chu S I 2011 Opt. Express 19 23857
[98] Li P C Zhou X X Wang G L Zhao Z X 2009 Phys. Rev. 80 053825
[99] Li P C Chu S I 2012 Phys. Rev. 86 013411